!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


module mod_psi_locmat
!General comments: computation of the (global) lumped mass matrix and
! of the local matrices: advection operator and forcing term.
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_base, only: &
   mod_base_initialized, &
   t_base

! use mod_grid, only: &
!   mod_grid_initialized, &
!   t_grid, grid,  &
!   t_e, affmap,   &
!   p_t_s, t_s
!
! use mod_testcases, only: &
!   mod_testcases_initialized, &
!   ndim,        &
!   coeff_diff,  &
!   coeff_adv,   &
!   coeff_re,    &
!   coeff_f,     &
!   coeff_neu,   &
!   coeff_rob
!
!!-----------------------------------------------------------------------
! 
! implicit none
!
!!-----------------------------------------------------------------------
!
!! Module interface
!
! public :: &
!   mod_psi_locmat_constructor, &
!   mod_psi_locmat_destructor,  &
!   mod_psi_locmat_initialized, &
!   lumped_mass, &
!   psi_locmat_nc
!   ! boundary conditions not required here (purely advective problem)
!
! private
!
!!-----------------------------------------------------------------------
!
!! Module types and parameters
!
!! Module variables
!
! ! public members
! logical, protected ::               &
!   mod_psi_locmat_initialized = .false.
! ! private members
! integer, parameter :: &
!   p1      = 1, &
!   nscheme = 2, &
!   psi     = 3, &
!   psi_pos = 4
! integer :: scheme_indx
! character(len=*), parameter :: &
!   this_mod_name = 'mod_psi_locmat'
!
!!-----------------------------------------------------------------------
!
!contains
!
!!-----------------------------------------------------------------------
!
! subroutine mod_psi_locmat_constructor(scheme)
!  character(len=*), intent(in) :: scheme
!
!  character(len=20+len_trim(scheme)) :: message
!  character(len=*), parameter :: &
!    this_sub_name = 'constructor'
!
!   !Consistency checks ---------------------------
!   if( (mod_messages_initialized.eqv..false.) .or. &
!       (mod_kinds_initialized.eqv..false.)    .or. &
!       (mod_base_initialized.eqv..false.)     .or. &
!       (mod_grid_initialized.eqv..false.)     .or. &
!       (mod_testcases_initialized.eqv..false.) ) then
!     call error(this_sub_name,this_mod_name, &
!                'Not all the required modules are initialized.')
!   endif
!   if(mod_psi_locmat_initialized.eqv..true.) then
!     call warning(this_sub_name,this_mod_name, &
!                  'Module is already initialized.')
!   endif
!   !----------------------------------------------
!
!   CaseName: select case(trim(scheme))
!    case("P1")
!     scheme_indx = p1
!    case("N-scheme")
!     scheme_indx = nscheme
!    case("PSI")
!     scheme_indx = psi
!    case("PSI+")
!     scheme_indx = psi_pos
!    case default
!     write(message,'(a,a,a)') 'Unknown scheme "',trim(scheme),'".'
!     call error(this_sub_name,this_mod_name,message)
!   end select CaseName
!
!   mod_psi_locmat_initialized = .true.
! end subroutine mod_psi_locmat_constructor
!
!!-----------------------------------------------------------------------
! 
! subroutine mod_psi_locmat_destructor()
!  character(len=*), parameter :: &
!    this_sub_name = 'destructor'
!   
!   !Consistency checks ---------------------------
!   if(mod_psi_locmat_initialized.eqv..false.) then
!     call error(this_sub_name,this_mod_name, &
!                'This module is not initialized.')
!   endif
!   !----------------------------------------------
!
!   mod_psi_locmat_initialized = .false.
! end subroutine mod_psi_locmat_destructor
!
!!-----------------------------------------------------------------------
! 
! pure subroutine lumped_mass(m)
!  real(wp), intent(out), allocatable :: m(:)
!
!  integer :: ie
!
!   allocate(m(grid%nv))
!   m = 0.0_wp
!   do ie=1,grid%ne
!     m(grid%e(ie)%iv) = m(grid%e(ie)%iv) + grid%e(ie)%a/3.0_wp
!   enddo
! 
! end subroutine lumped_mass
! 
!!-----------------------------------------------------------------------
! 
! pure subroutine psi_locmat_nc( aak, fk, base, e , u )
! 
!  real(wp), intent(out) :: aak(:,:), fk(:)
!  real(wp), intent(in) :: u(:)
!  type(t_base), intent(in) :: base
!  type(t_e), intent(in) :: e
! 
!  logical, parameter :: & ! masks to classify the element
!    mask1t(3) = (/ .true. , .false. , .false. /), & ! one target el
!    mask2t(3) = .not.mask1t                         ! two target el
!  logical :: kk_gt0(3), onetarget
!  integer :: l, i, ii(3)
!  real(wp) :: a(ndim,base%m), f(base%m), xyg(ndim,base%m), &
!    abar(ndim), kgphi(ndim,3), kk(3), beta(3)
!
!   ! Gaussian nodes in physical space
!   xyg = affmap(e,base%xig)
!
!   ! problem coefficients
!   a = coeff_adv(xyg)
!   f = coeff_f(xyg)
!
!   ! average a
!   abar = 0.0_wp
!   do l=1,base%m
!     abar = abar + base%wg(l)*a(:,l)
!   enddo
!   abar = 2.0_wp*abar
!
!   ! Compute the coefficients
!   !   kk_i = |K| a \cdot \nabla\phi_i
!   ! using the fact that
!   !   |K|\nabla\phi_i = 1/2 |e_i|n^{int}_i
!   ! and n^{int} is the iterior normal.
!   do i=1,3
!     kgphi(:,i) = -0.5_wp * e%s(i)%p%l * e%n(:,i)
!     kk(i) = dot_product( abar , kgphi(:,i) )
!   enddo
!
!   ! Classify the element
!   classify: do i=1,3
!     ii(1) = i; ii(2:3) = mod(ii(1)+(/0,1/),3)+1
!     kk_gt0 = kk(ii) .gt. 0.0_wp ! test where kk.gt.0
!
!     if( all(kk_gt0.eqv.mask1t) ) then
!       ! found a one target element
!       onetarget = .true.
!       ! ii stores the indexes as   ii = (/ out , in , in /)
!       exit classify
!     elseif( all(kk_gt0.eqv.mask2t) ) then
!       ! found a two target element
!       onetarget = .false.
!       ! ii stores the indexes as   ii = (/ in , out , out /)
!       exit classify
!     endif
!   enddo classify
!
!   CaseName: select case(scheme_indx)
!    case(p1)      ! P1 non-monotone
!     call p1_amat(aak,kk)
!    case(nscheme) ! N-scheme
!     call n_scheme_amat(aak,onetarget,ii,kk)
!    case(psi)     ! PSI
!     call psi_amat(aak,beta,onetarget,ii,kk,u)
!    case(psi_pos) ! PSI with positive matrix
!     call psi_amat_pos(aak,beta,onetarget,ii,kk,u)
!   end select CaseName
!
!   ! the betas of the N-scheme are not bounded, so we can't use them
!   ! to off-center the source term
!   fk  = 0.0_wp
!   do l=1,base%m
!     do i=1,3
!       fk(i) = fk(i) + base%wg(l) * f(l) * base%p(i,l)
!     enddo
!   enddo
!   fk  = 2.0_wp*e%a*fk
! 
! end subroutine psi_locmat_nc
! 
!!-----------------------------------------------------------------------
! 
! pure subroutine p1_amat(aak,kk)
!  real(wp), intent(out) :: aak(3,3)
!  real(wp), intent(in) :: kk(3)
!
!  integer :: i
!
!   do i=1,3
!     aak(i,:) = kk/3.0_wp
!   enddo
! 
! end subroutine p1_amat
! 
!!-----------------------------------------------------------------------
! 
! pure subroutine n_scheme_amat(aak,onetarget,ii,kk)
!  real(wp), intent(out) :: aak(3,3)
!  logical, intent(in) :: onetarget
!  integer, intent(in) :: ii(3)
!  real(wp), intent(in) :: kk(3)
!
!   onetwo_targets: if(onetarget) then
!     aak(ii(1),:) = kk
!     aak(ii(2),:) = 0.0_wp
!     aak(ii(3),:) = 0.0_wp
!   else
!     aak(ii(1),: ) = (/   0.0_wp   ,  0.0_wp   ,  0.0_wp  /)
!     aak(ii(2),ii) = (/ -kk(ii(2)) , kk(ii(2)) ,  0.0_wp  /)
!     aak(ii(3),ii) = (/ -kk(ii(3)) ,  0.0_wp   , kk(ii(3))/)
!   endif onetwo_targets
! 
! end subroutine n_scheme_amat
! 
!!-----------------------------------------------------------------------
!
! pure subroutine psi_amat(aak,beta,onetarget,ii,kk,u)
! ! Here we also compute the betas, which can be useful for the
! ! off-centering of the source term
!  real(wp), intent(out) :: aak(3,3), beta(3)
!  logical, intent(in) :: onetarget
!  integer, intent(in) :: ii(3)
!  real(wp), intent(in) :: kk(3), u(3)
! 
!  integer :: i,j, iminmax(2)
!  real(wp) :: sgrad(2), res, lambda
!
!   onetwo_targets: if(onetarget) then
!     call n_scheme_amat(aak,.true.,ii,kk)
!     beta(ii) = (/ 1.0_wp , 0.0_wp , 0.0_wp /)
!   else
!     beta(ii(1)) = 0.0_wp
!     sgrad = u(ii(2:3)) - u(ii(1)) ! side gradients
!     lambda = product(sgrad) ! determinant
!     res = sum( kk(ii(2:3))*sgrad ) ! residual
!     if(lambda.ge.0.0_wp) then
!       call n_scheme_amat(aak,.false.,ii,kk)
!       if(res.ne.0.0_wp) then
!         beta(ii(2:3)) = kk(ii(2:3))*(sgrad/res)
!       else
!         beta(ii(2:3)) = 0.5_wp ! not uniquely defined
!       endif
!     else
!       ! find the maximum among ii(2) and ii(3)
!       if(u(ii(2)).le.u(ii(3))) then
!         iminmax = (/ 2, 3 /)
!       else
!         iminmax = (/ 3, 2 /)
!       endif
!       ! now assign the residual
!       if(res>0.0_wp) then ! residual tends to diminsh the solution
!         beta(ii(iminmax)) = (/ 0.0_wp , 1.0_wp /)
!       else ! residual tends to increase the solution
!         beta(ii(iminmax)) = (/ 1.0_wp , 0.0_wp /)
!       endif
!       do i=1,3
!         do j=1,3
!           aak(i,j) = beta(i)*kk(j)
!         enddo
!       enddo
!
!     endif
!   endif onetwo_targets
! 
! end subroutine psi_amat
! 
!!-----------------------------------------------------------------------
!
! pure subroutine psi_amat_pos(aak,beta,onetarget,ii,kk,u)
! ! Analogous to psi_amat, but computes and equivalent matrix which is
! ! positive. The reason for introducing a separate function is that in
! ! some special cases the positive matrix of the PSI scheme is not
! ! well defined, and so we have to resort to the N-scheme. In these
! ! cases, however, the error caused by N-scheme is very small.
!  real(wp), intent(out) :: aak(3,3), beta(3)
!  logical, intent(in) :: onetarget
!  integer, intent(in) :: ii(3)
!  real(wp), intent(in) :: kk(3), u(3)
! 
!  integer :: i,j, iminmax(2), i0, i1
!  real(wp) :: sgrad(2), res, lambda
!  real(wp) :: delta, coeff, kkmod(3)
!
!   onetwo_targets: if(onetarget) then
!     call n_scheme_amat(aak,.true.,ii,kk)
!     beta(ii) = (/ 1.0_wp , 0.0_wp , 0.0_wp /)
!   else
!     beta(ii(1)) = 0.0_wp
!     sgrad = u(ii(2:3)) - u(ii(1)) ! side gradients
!     lambda = product(sgrad) ! determinant
!     res = sum( kk(ii(2:3))*sgrad ) ! residual
!     if(lambda.ge.0.0_wp) then
!       call n_scheme_amat(aak,.false.,ii,kk)
!       if(res.ne.0.0_wp) then
!         beta(ii(2:3)) = kk(ii(2:3))*(sgrad/res)
!       else
!         beta(ii(2:3)) = 0.5_wp ! not uniquely defined
!       endif
!     else
!       ! find the maximum among ii(2) and ii(3)
!       if(u(ii(2)).le.u(ii(3))) then
!         iminmax = (/ 2, 3 /)
!       else
!         iminmax = (/ 3, 2 /)
!       endif
!       ! now assign the residual
!       if(res>0.0_wp) then ! residual tends to diminsh the solution
!         i0 = 1; i1 = 2
!       else ! residual tends to increase the solution
!         i0 = 2; i1 = 1
!       endif
!       ! At this point we could already build the matrix, and the
!       ! method would be monotone. The matrix however would not be
!       ! monotone, so that it could not be used to get a steady state
!       ! solution with an implicit time-stepping. It's possible
!       ! however to modify the matrix obtaining \tilde{A} such that
!       !  \tilde{A} u = A u
!       ! for the present value u (but not for a generic vector) and
!       !  \tilde{A} is positive
!       ! \tilde{A} has to be computed adding a correction to u which
!       ! is orthogonal to Grad(u), so that the residual does not
!       ! change.
!       beta(ii(iminmax((/i0,i1/)))) = (/ 0.0_wp , 1.0_wp /)
!       delta = -sgrad(iminmax(i1)-1)
!       if(abs(delta).lt.10.0_wp*kk(ii(iminmax(i0)))*epsilon(delta)) then
!         call n_scheme_amat(aak,.false.,ii,kk)
!         return
!       endif
!       coeff = kk(ii(iminmax(i0)))/delta
!       kkmod(ii(1)) = kk(ii(1)) - &
!                      coeff*(u(ii(iminmax(i1)))-u(ii(iminmax(i0))))
!       kkmod(ii(iminmax(i1))) = kk(ii(iminmax(i1))) - &
!                      coeff*(u(ii(iminmax(i0)))-u(ii(1)))
!       kkmod(ii(iminmax(i0))) = 0.0_wp
!       do i=1,3
!         do j=1,3
!           aak(i,j) = beta(i)*kkmod(j)
!         enddo
!       enddo
!
!     endif
!   endif onetwo_targets
! 
! end subroutine psi_amat_pos
! 
!!-----------------------------------------------------------------------

end module mod_psi_locmat

